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An Introduction to Using Software Tools 
for Automatic Differentiation 



by 

Uwe Naumann and Andrea Walther 

Abstract 

We give a gentle introduction to using various software tools for automatic differentiation 
(AD). Ready-to-use examples are discussed, and links to further information are presented. Our 
target audience includes all those who are looking for a straightforward way to get started using 
the available AD technology. The document is dynamic in the sense that its content will be 
updated as the AD software evolves. 

1 Intention and References 

This document explains how to use the following software tools for automatic differentiation (AD) 
to generate first-order derivative code for a small example problem. 

• ADIFOR 2.0 Revision D (http://www.mcs.anl.gov/adifor) 

• ADIC 1.1 (http://www.mcs.anl.gov/adic) 

• ADOL-C 1.8 (http://www.math.tu-dresden.de/wir/project/adolc) 

• TAPENADE 1.0-alpha j |http : //www-sop . inria . f r/tropics| | 

The document does not describe the AD algorithms used by these software tools. Various pub- 
lications cover the theoretical foundations of the algorithms referenced in our discussion. For an 
introduction and more advanced issues we point the reader to Three workshops have been 
dedicated to AD, and the proceedings |SJ E] contain numerous references to other AD-related 
papers that have appeared in scientific journals and proceedings of international conferences. A 
large variety of successful applications of AD software to real- world applications are discussed there 
as well. 

2 Getting Started 

For information on how to obtain, install, and use the AD tools discussed in this document, visit 
the Web sites listed in Section 1. The examples in this document are discussed in the context of the 
computer environment found on a laptop computer, specifically an Intel Pentium system running 
Mandrake Linux Version 8.1. ADIFOR, ADIC, and TAPENADE are installed in the directory 
/home/uwe/ADTODLS. The ADOL-C library and include files are assumed to be in the same directory 
as the source files. 

ADIFOR, ADIC, and TAPENADE are available on the Web and can be accessed online by using 
the corresponding Web servers. Visit 

• http : //www.mcs . anl .gov/ autodif f /adif or server for ADIFOR, 

• |http: //www.mcs .an l .gov/autodif f /adicserver|for ADIC, and 

• http ://tapenade. inria. fr:8080/tapenade/index. jsp for TAPENADE. 
All source files presented in this document can be downloaded from 

http : //www-unix .mcs . anl .gov/ "naumann/ ad_tools .html 
The makefiles have to be adapted to the user's environment. 
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3 Explosion Equation 

We have selected a variant of the Bratu problem pp to illustrate the use of source transformation AD 
tools. The code in Appendix lAl models the thermal explosion of solid fuels, which can be described 
by the system of differential equations 

x"{t) + s • e = 0, 

where r £ (—1, 1) and x(— 1) = x(l) = 0. The problem has been discretized by using step size h as 

Fi = Xi-i - 2 Xi + x i+1 + /i 2 [/i_i + 10/, + /i+i]/12 

for i = 1, . . . , 10000, with xo — £10001 — and fi = s exp(xi/(l+tXi)). Of interest are the derivatives 
of the component functions Fi with respect to the current state Xi as well as the parameters s and 
t. 

We consider the following problems: 

1. ADIFOR 

(a) Use the forward vector mode to compute the Jacobian of F with respect to x and the 
parameters s and t. In our implementation, F is represented by the variable f , x is repre- 
sented by x, and s and t are combined into the parameter vector prm. All three program 
variables are declared as arrays of double-precision floating-point numbers, namely, 

• double precision x(7) , prm(2) , f(dim) in Fortran and 

• double x[7], prm [2] , f [7] in C. 

(See Section 0) 

(b) Use Curtis-Powell- Rcid 8 seeding to compute the compressed Jacobian in forward vector 
mode (see Section 13.1. 1(1 . 

(c) Use the SparsLinC library to compute the sparse Jacobian by sparse forward mode (see 
Section 13. 1.2(1 . 

2. ADOL-C 

(a) Use one of the easy-to-use drivers to compute the Jacobian of F with respect to x and 
the parameters s and t. 

(b) Compute the sparsity structure of the Jacobian. 

(c) Use the low-level routines for forward and reverse mode to compute the Jacobian of F 
with respect to x and the parameters s and t. 

(See Section El) 

3. ADIC 

(a) Compute the Jacobian of F with respect to x in forward vector mode (see Section l3~2l . 

4. TAPENADE 

(a) Compute the Jacobian of F with respect to both x and the parameters s and t in reverse 
mode (see Section l3~lj) . 



2 



3.1 ADIFOR 2.0 

ADIFOR generates a differentiated version of the subroutine expl (...) with the following header. 

subroutine g_expl(g_p_, dim, parmax, x, g_x, ldg_x, prm, g_prm, 
+ ldg_prm, f, g_f, ldg_f) 
integer dim, parmax 

double precision x(dim), prm(parmax) 
double precision f (dim) 

integer g_pmax_ 

parameter (g_pmax_ = 9) 

integer g_p_ , ldg_x, ldg_prm, ldg_f 

double precision g_x(ldg_x, dim), g_prm(ldg_prm, parmax), 
+ g_f(ldg_f, dim) 

The parameter g_pmax_ is the maximal number of directional derivatives that can be computed. Its 
value is set by using AD_PMAX in explosion, adf. The parameter is required because of the lack of 
dynamic memory allocation in Fortran 77. The derivative components of all scalar program variables 
are allocated as vectors of length g_pmax_. For example, in order to compute the Jacobian matrix 
of f with respect to x, the value of g_pmax_ must be at least equal to 7, so that the derivative 
components of x can accommodate the identity in JR 7 . 

The actual number of directional derivatives g_p_ must be less than or equal to g_pmax_. The 
parameter g_p_ determines the upper bound for the loops that compute the values of the derivative 
components. 

To compute the Jacobian of f with respect to both x and prm, we set g_pmax_ to be greater than 
or equal to the sum of the numbers of elements in both vectors, that is, 9 = 7 + 2 = dim + parmax. 
Consequently, g_x contains the first dim columns of the seed matrix and g_prm its last parmax 
columns. The argument g_f is the transpose of the (dim x ldgJ:)- Jacobian, where ldgj must be 
initialized to 9. All this is done in the driver routine, which is shown in Appendix IB. 3. 21 

The "dense:" section of the makefile in Appendix IB. 21 builds an executable that produces the 
following result. 
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The output has been formatted for better readability. It shows the full 7x9 Jacobian evaluated at 
the argument: 

x(l) = 1.72 

x(2) = 3.45 

x(3) = 4.16 

x(4) = 4.87 

x(5) = 4.16 

x(6) = 3.45 

x(7) = 1.72 
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prm(l) = 1.3 
prm(2) = 0.245828 



3.1.1 Compressed Jacobian Seeding 



The full Jacobian computed in the preceding section is sparse, and its sparsity pattern can be 
visualized as follows. 



* 

*. 

* 



Here, * stands for a nonzero entry, and blanks represent structural zero entries in the Jacobian. In 
other words, no dependence exists between the corresponding dependent and independent variables. 

Curtis-Powell-Rcid (CPR) |S] seeding is based on the idea that certain columns of the Jacobian 
can be merged to share storage. For example, column 1 and column 7 could share one column, 
thus resulting in a compressed version of the Jacobian. This implies that the sparsity pattern must 
be known in advance in order to exploit matrix compression techniques. Recall that the derivative 
code generated by ADIFOR always loops over the derivative components of all active variables.* 
Many of them are equal to zero, leading to predictably trivial multiplications that one would like to 
avoid. Therefore, instead of computing the Jacobian as a Jacobian times identity matrix product, 
one could try to compute a compressed Jacobian using a seed matrix with fewer columns than the 
identity. Since the number of independent variables is often very large, the size of the seed matrix 
can become much smaller, leading to a decreased complexity of the forward vector mode Jacobian 
computation. 

In CPR seeding, one considers the column incidence graph of the Jacobian to try to determine 
a minimal vertex coloring. Whenever two vertices share the same color, the corresponding columns 
can be stored in the same column of the compressed Jacobian. In our example, the column incidence 
graph has the following structure. 





y 




00 









+An active variable w can be characterized as follows. At some point in the program (1) the value of w depends on 
the value of some independent variables and (2) there is some dependent variable whose value depends on w. Variables 
that are not active are referred to as passive. 
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Unfortunately, since the vertex coloring problem is known to be NP-complete in general, the use of 
heuristics is essential. The coloring in the example graph has been found "by inspection." Different 
colors are represented by different vertex shapes. 

The number v of different colors used determines the number of columns in the CPR seed 
matrix. Its rows are Cartesian basis vectors in M v . Whenever two vertices share the same color, 
the corresponding rows in the seed matrix contain the same Cartesian basis vector. This leads to 
the following seed matrix for our example. 

1 
10 

10 

1 
10 

10 

1 
1 
1 

The "compressed:" section of the makefile in Appendix IB. 21 builds an executable that produces 
output similar to the following. 
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It shows the compressed 7x5 Jacobian evaluated at the current argument. The reconstruction of 
the original Jacobian is a simple substitution process as described in Chapter 7]. 

3.1.2 Sparse Forward Mode SparsLinC 

We compute the Jacobian of f with respect to x. The sparsity pattern of the Jacobian need not be 
known a priori. Derivative components are sparse vectors of (index, value) pairs. The computational 
overhead of sparse vector arithmetic results from the index calculations. Structural sparsity of the 
extended Jacobian 9 , Chapter 2] is exploited by avoiding trivial multiplications by zero. The decision 
about when to apply runtime sparsity methods depends on the problem structure. 
The use of SparsLinC with our example can be described by the following steps: 

1. Use the file explosion. cmp without change. 

2. Add the entry AD_FLAVDR=sparse to explosion . adf. 

3. Write driver program as shown in Appendix IB. 51 

4. Generate g_explosion.f by calling Adif or2. 1 AD_SCRIPT=explosion . adf , and copy it from 
output Jiles to the current directory. 

5. Compile g_explosion . f and the driver program. 

6. Link with ReqADIntrinsics-Limix86 . o, libADIntrinsics-Linux86 . a, and 
libSparsLinC-Linux86 . a. 
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7. Run the executable. 

Let us have a closer look at the driver program. The differentiated subroutine 

g_expl (dim , parmax , x , g_x , prm , f , g J ) 

is generated in g_explosion . f , where g_x and f_x are integer arrays of dimension dim containing 
pointers to the corresponding sparse derivative objects. Seeding g_x is performed in the driver by 

do 10 i=l,dim 
g_x(i)=0 

CALL DSPSD(g_x(i) ,i,l.d0,l) 
10 continue 

Notice that g_x(i) must be initialized properly before calling DSPSD. The call initializes the sparse 
derivative object pointed at by g_x(i) as (i.l.dO). Consequently, after executing this loop, g_x 
contains the sparse identity in J£ dlm . 

The sparse derivative components are extracted by 

do 30 i=l,dim 
g_f (i)=0 

CALL DSPXSQ (indvec , valvec ,dim,g_f (i) ,outlen, info) 
30 continue 

where indvec is the index vector and valvec the corresponding value vector making up the sparse 
representation of the derivative object pointed at by g_f (i) . The parameter outlen is the number of 
nonzero entries in the 1th row of the Jacobian. Consequently, the outlen first entries of indvec and 
valvec define the nonzero entries in the 1th row of the Jacobian as index value pairs (indvec(j) , 
valvec ( j ) ). 

The "sparse:" section of the makefile in Appendix IB . 21 builds an executable that produces out- 
put similar to the following. 
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The sparse Jacobian is given by sparse row vectors whose nonzero entries are represented by (index, 
value) pairs. The reconstruction of the full Jacobian is straightforward. Refer to |3] for further 
information on how to use SparsLinC. 

3.2 ADIC 1.1 

All files involved in using ADIC with the C version of our example problem are shown in Appendix ICl 
The Jacobian of F with respect to x is computed in forward vector mode. ADIC uses the derived 
data type DERIV_TYPE to associate derivative components with active variables. The data type 
InactiveDouble is used to indicate that some floating-point variable is not active. The call of 
ad_AD_Init (dim) causes the derivative code to use only the first dim elements of the derivative 
components. The latter are vectors of length GRAD_MAX (see Appendix EJ- The vector x is declared 
to contain the dim independent variables by calling 
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ad_AD_SetIndepArray(x,dim) ; 
ad_AD_SetIndepDone() ; 

The function value components of DERIV_TYPE variables are accessed by means of the macro DERIV_val. 

ADIC generates a differentiated version of the subroutine explosion and names it ad_explosion. 
The argument list remains unchanged because all AD-related information is encapsulated in the new 
data type DERIV_TYPE. If the function returns a double, however, the differentiated function becomes 
a procedure, and the returned value becomes the first argument. The derivative components of scalar 
variables of this type can be extracted by calling the routine ad_AD_ExtractGrad. The actual function 
value can be extracted with ad_AD_ExtractVal. In the example we call ad_AD_ExtractGrad( jac ,F [i] ) 
inside a loop over the dim elements of the vector of dependent variables F. The auxiliary variable 
jac is declared as a passive vector of size dim and contains the 1th row of the Jacobian. 

A makefile is provided to build the executable explosion, ad in order to compute the first seven 
columns of the Jacobian shown in Section 13.11 

3.3 ADOL-C 

The AD-tool ADOL-C is based on operator overloading. Using this technique, one can log for each 
operation during the program execution the operator and the variables that are involved. Hence, one 
obtains a new internal representation of the function evaluation. Based on the generated execution 
log, ADOL-C computes the desired derivatives. 

To apply ADOL-C for derivative calculations, one first has to modify the evaluation program 
to record the internal representation called tape. This modification starts with including header- 
file^) that introduce the new data types and functions. Here, the easiest way is to simply include 
adolc.h. Second, one has to define the part of the program for which one wants to compute the 
derivatives. From now on, this part is called the active section. The statement trace_on(tag,fceep); 
determines the beginning of the active section, the statement trace_off(/iZe); the end of the active 
section. The parameter tag identifies the function to be differentiated. Hence, several function 
representations can be kept at the same time. Because of the shortness of this introduction, the 
optional parameters keep and file are not explained here but are explained in the documentation 0]. 
Third, one has to change the types of the independents to adoubles and mark them as independents 
using the overloaded operator <<=. Similarly, one must mark the dependents using the overloaded 
operator >>=. Finally, all variables that lie on the way from the independents to the dependents 
must be declared as adoubles. This step includes the generation of a new function explosion_ad, 
which contains the same source code as before but in the interface the double-variables have to be 
changed to adouble- variables. The required modification of the original source code is illustrated by 
Appendix IdI 

3.3.1 Jacobian Calculation with the Easy-to-Use Drivers 

After the tapes are generated during the execution of the active section, the required derivative 
objects can be calculated. For that purpose ADOL-C provides a variety of easy-to-use drivers: 
gradient(..), jacobian(..), hessian(..), vec_jac(..) computing the product vector times Jacobian, 
jac_vec(..) computing the product Jacobian times vector, and so on. Moreover, ADOL-C supplies 
routines evaluating the Taylor coefficient vectors and their Jacobians with respect to the current 
state vector of solution curves defined by ordinary differential equations. Furthermore, there are 
drivers for derivative tensors and for the differentiation of implicit and inverse functions. 

If one wants to compute the Jacobian of F with respect to x and the parameters s and t, one 
has to declare a variable storing the derivative information 

double **J = myalloc(dim,dim+parmax); 
where myalloc is provided by ADOL-C to allocate a two-dimensional array. Then, the statement 
jacobian (tag,dim,dim+parmax,v, J); 
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causes the computation of the Jacobian of the function representation contained in the tape with 
the number tag at the point v. For a consistency check, the second and third parameter determine 
the number of dependents and independents, respectively. Hence, only the two statements given 
above have to be inserted after trace_off() in order to compute the Jacobian of F. 

3.3.2 Calculation of Sparsity Pattern 

ADOL-C provides for the computation of sparsity patterns the driver 

jac_pat(tag,dim 1 dim+parmax,v,rb,cb,Jsp, option); 
If one sets rb and cb to NULL, jac_pat computes the sparsity structure of the complete Jacobian at the 
point where the tape was generated and stores it in the unsigned int-array Jsp. The corresponding 
statements that have to be added to the source code are shown in Appendix D. If a certain block 
structure of the Jacobian is known, the unsigned int-vectors rb and cb can be applied to describe a 
compressed form of the independent and dependent variables. 

3.3.3 Forward and Reverse Mode Using Low-Level Routines 

Arbitrarily high-order derivatives can be calculated by using the low-level functions of ADOL-C 
for the forward and reverse mode of AD. These routines are explained in detail in a short reference 
available from the ADOL-C Web page. In this article, only the computation of first-order derivatives 
is sketched. If one wants the full Jacobian, it is preferable to use vector modes of AD. For that 
purpose, ADOL-C provides the drivers fov_forward(...) and fov_reverse (...). Here, fov stands for 
first-order vector. The other acronyms of the low-level routines have a corresponding meaning. The 
statement 

fov_forward(tag,dim,dim+parmax,p,v,X,Fp,Y); 

computes the derivative object Y = F'(v)X for X E j^im+parmaxxp Thc statement 

fov_reverse(tag,dim,dim+parmax,q,U,Z); 
computes the derivative object Z = l)F'(v) for U E M qxd,m . To prepare this reverse sweep, one 
has to call an appropriate forward routines, the choice of which is described in detail by the short 
reference mentioned above. The code segments required for computing the full Jacobian with the 
low- level routines are contained in Appendix D. 

3.4 TAPENADE 

Our last example uses the alpha version of TAPENADE 1.0 to illustrate the use of reverse-mode 
AD for computing the Jacobian of f with respect to x and prm. All files needed to run the example 
are described in Appendix El The makefile in Appendix IE. II can be used to generate the executable 
explosion . ad. Running the latter results in the computation of the full Jacobian as in Section 13. II 
but this time using the reverse mode of AD. 

The following command-line parameters are used to call tapenade: 

• "-head" : The name of the head routine 

• "-vars" : The names of the independent variables 

• "-cl": Switch indicating that reverse, or cotangent linear, mode is used 

It generates derivative code in explcl . f , which must be compiled together with the driver program 
shown in Appendix IE. 21 and linked with the routines for storing and restoring the values of the tape 
as described in Chapter 2]. 

The current version of TAPENADE does not provide the reverse vector mode. Hence, thc 
Jacobian must be computed as a sequence of Jacobian transposed times vector products. This 
computation is done in the driver program by successively initializing the derivative components 
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g_f of f as the Cartesian basis vectors in M dlm . Thus, the Jacobian is accumulated row by row 
at a computational complexity proportional to the number of dependent variables. This feature 
becomes particularly interesting in the case of large, single gradients. Refer to [5] for further details 
on reverse-mode AD. 
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A Original Code 



subrout ine bratu (dim , parmax , x , prm , F) 

integer dim , parmax 
C independent variables 

double precision x(dim) , prm(parmax) 
C dependent variables 

double precision F(dim) 

C 

integer i 

double precision h 
h = 2.0/(dim+l) 

F(l) = -2*x(l)+h*h*prm(l)/12.0*(l+10*exp(x(l)/(1.0+prm(2)*x(l)))) 
F(2) = x(l)+h*h*prm(l)/12.0*exp(x(l)/(1.0+prm(2)*x(l))) 

do 1 i=2,dim-l 

F(i-l) = F(i-l)+x(i)+h*h*prm(l)/12.0*exp(x(i)/(1.0+prm(2)*x(i))) 
F(i) = F(i)-2*x(i)+h*h*prm(l)/1.2*exp(x(i)/(1.0+prm(2)*x(i))) 
F(i+1) = x(i)+h*h*prm(l)/12.0*exp(x(i)/(1.0+prm(2)*x(i))) 
1 continue 

F(dim-l) = F(dim-l)+x(dim)+h*h*prm(l)/12.0*exp(x(dim)/(1.0 

* +prm(2)*x(dim))) 
F(dim) = F(dim)-2*x(dim) 

F(dim) = F(dim)+h*h*prm(l)/12.0*(l+10*exp(x(dim)/(1.0 

* +prm(2)*x(dim)))) 
end 

B ADIFOR 2.0 

The following files can be downloaded from 

: //www-unix .mcs . anl . gov/~naumann/ ad_tools .html 

explosion. adf 
explosion. cmp 

explosion. driver. compressed. f 
explosion . driver . f 
explosion. driver . sparse . f 
explosion. f 
explosion . sparse . adf 
makefile 

B.l explosion. cmp 

The composition file lists the names of all files containing subroutines subject to differentiation, 
our simple example there is just explosion. f. 
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B.2 makefile 



The makefile can be used for computing the full, compressed, and sparse Jacobians using ADI- 
FOR 2.0. It also ensures a proper cleanup of all files that are generated automatically during this 
process. 

AD_LIB=/home/uwe/ ADT00LS/ ADIF0R2/ ADIFDR2 . OD . lib 
dense : 

Adif or2 . 1 AD_SCRIPT=explosion . adf 

cp output_f iles/g_explosion.f . 

g77 -g -c g_explosion . f explosion . driver . f 

g77 -g -o explosion. ad. dense -L$ (AD_LIB) /lib *.o \ 

$ (AD_LIB) /lib/ReqADIntrinsics-Linux86 . o \ 

-lADIntrinsics-Linux86 

compressed: 

Adif or2 . 1 AD_SCRIPT=explosion . adf 
cp output_f iles/g_explosion.f . 

g77 -g -c g_explosion . f explosion . driver . compressed. f 
g77 -g -o explosion . ad . compressed -L$(AD_LIB)/lib *.o \ 
$ (AD_LIB) /lib/ReqADIntrinsics-Linux86 . o \ 
-lADIntrinsics-Linux86 

sparse : 

Adif or2 . 1 AD_SCRIPT=explosion . sparse . adf 
cp output_f iles/g_explosion.f . 

g77 -g -c g_explosion . f explosion . driver . sparse . f 
g77 -g -o explosion . ad . sparse -L$ (AD_LIB) /lib *.o \ 
$ (AD_LIB) /lib/ReqADIntrinsics-Linux86 . o \ 
-lADIntrinsics-Linux86 \ 
-lSparsLinC-Linux86 

clean: 

rm -fr output_f iles 
rm -fr AD_cache 
rm * . o 
rm g_* 

rm explosion. ad. * 
rm *\~ 

B.3 Jacobian 
B.3.1 explosion. adf 

The Jacobian of the output variable f with respect to the two input variables x and prm of the 
top-level routine expl is computed. The parameter AD_PMAX must be set to the total number of 
independent variables, that is, 9 § . 

AD_PR0G = explosion. cmp 

§This is because the Jacobian is computed by forward vector mode with a seed matrix that is equal to the identity 
in Ft 9 . 
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AD_TOP = expl 
AD_PMAX = 9 
AD_IVARS= x prm 
AD_DVARS= f 

B.3.2 explosion. driver. f 

program main 
implicit none 

C 

C Example: Explosion Equation 
C Driver for computing Jacobian 
C 

integer dim, parmax, n 

parameter (dim=7, parmax=2, n=9) 

integer i,j 
C independent variables 

double precision x(dim) , prm(parmax) 
C derivative components of independent variables 

double precision g_x(n,dim) 

double precision g_prm(n, parmax) 
C dependent variables 

double precision f (dim) 
C derivative components of dependent variables 

double precision g_f (n,dim) 

C Initialization of input variables 
x(l) = 1.72 
x(2) = 3.45 
x(3) = 4.16 
x(4) = 4.87 
x(5) = 4.16 
x(6) = 3.45 
x(7) = 1.72 
prm(l) = 1.3 
prm (2) = 0.245828 

C Seeding (identity) 
do 20 i=l,n 

do 10 j=l, parmax 

if (i.eq.j+dim) then 

g_prm(i, j)=1.0 
else 

g_prm(i, j)=0.0 
endif 
10 continue 
20 continue 

do 40 i=l,n 
do 30 j=l,dim 

if (i.eq.j) then 
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g_x(i,j)=1.0 
else 

g_x(i,j)=0.0 
endif 
30 continue 
40 continue 

C call differentiated subroutine 

call g_expl (n , dim , parmax , x , g_x , n , prm , g_prm , n , f , g_f , n) ; 

C print Jacobian 
do 60 i=l,dim 
do 50 j=l, n 

print*, "f'(", i, j, ")=",g_f (j,i) 

50 continue 
60 continue 

end 

B.4 Compressed Jacobian 
B.4.1 explosion. adf 

This is the same as in the dense case. 

B.4. 2 explosion. driver. compressed. f 

program main 
implicit none 

C 

C Example: Explosion Equation 

C Driver for computing compressed Jacobian 

C 

integer dim, parmax, n 

parameter (dim=7, parmax=2, n=5) 

integer i,j 
C independent variables 

double precision x(dim), prm(parmax) 
C derivative components of independent variables 

double precision g_x(n,dim) 

double precision g_prm(n, parmax) 
C dependent variables 

double precision f (dim) 
C derivative components of dependent variables 

double precision g_f (n.dim) 

C Initialization of input variables 
x(l) = 1.72 
x(2) = 3.45 
x(3) = 4.16 
x(4) = 4.87 
x(5) = 4.16 
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x(6) = 3.45 
x(7) = 1.72 
prm(l) = 1.3 
prm(2) = 0.245828 

C Seeding (CPR) 
do 20 i=l,n 

do 10 j=l,parmax 
g_prm(i, j)=0.0 
10 continue 
20 continue 

g_prm(n-l,l)=1.0 
g_prm(n,2)=1.0 

do 40 i=l,n 
do 30 j=l,dim 
g_x(i,j)=0.0 
30 continue 
40 continue 

g_x(l,l)=1.0 
g_x(l,4)=1.0 
g_x(l,7)=1.0 
g_x(2,2)=1.0 
g_x(2,5)=1.0 
g_x(3,3)=1.0 
g_x(3,6)=1.0 

C call differentiated subroutine 

call g_expl (n , dim , parmax , x , g_x , n , prm , g_prm , n , f , g_f , n) ; 

C print compressed Jacobian 
do 60 i=l,dim 
do 50 j=l,n 

print*, "f'(", i, j, ")=",g_f (j,i) 

50 continue 
60 continue 

end 

B.5 Sparse Jacobian 
B.5.1 explosion. sparse. adf 

In order to make ADIFOR generate derivative code that can use the sparse forward mode provided 
by SparsLinC, the parameter AD_FLAV0R must be set to sparse. 

AD_PR0G = explosion. cmp 
AD_FLAV0R = sparse 
AD_T0P = expl 
AD_PMAX = 9 
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AD_IVARS= x 
AD_DVARS= f 



B.5.2 explosion, driver, sparse, f 

program main 
implicit none 

C 

C Example: Explosion Equation 

C Driver for computing sparse Jacobian using SparsLinC 
C 

integer dim, parmax, n 

parameter (dim=7, parmax=2, n=9) 

integer i,j 
C independent variables 

double precision x(dim) 
C passive inputs 

double precision prm(parmax) 
C pointers to sparse derivative components of 
C independent variables 

integer g_x(dim) 
C dependent variables 

double precision f (dim) 
C pointers to sparse derivative components of 
C dependent variables 

integer g_f (dim) 

C (index, value) pairs 

integer indexes (dim) 
double precision values (dim) 

C values used in extraction routine 
integer outlen, info 



C Initialization of input variables 
x(l) = 1.72 
x(2) = 3.45 
x(3) = 4.16 
x(4) = 4.87 
x(5) = 4.16 
x(6) = 3.45 
x(7) = 1.72 
prm(l) = 1.3 
prm(2) = 0.245828 

C Initialization of sparse data structures 
call XSPINI 
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c 



Seeding (sparse identity) 



do 10 i=l,dim 
g_x(i)=0 
g_f (i)=0 

call DSPSD(g_x(i) ,i,l.d0,l) 
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continue 



C 



Call differentiated subroutine 



call g_expl (dim , parmax , x , g_x , prm , f , g_f ) 



C 



Extract derivative components 



do 30 i=l,dim 

call DSPXSQ (indexes .values ,dim,g_f (i) , outlen, info) 
if (info.eq.O) then 
do 20 j=l, outlen 

print*, "indexes (", i, j, ")=", indexes(j) 

print*, "values (", i, j, ")=", values(j) 

20 continue 

endif 
30 continue 

end 



explosion. c 
explosion . driver . c 
explosion. init 
makefile 

C.l makefile 

AD_INC = -I$(ADIC)/include -I. 
AD_LIB = -L$(ADIC)/lib/$(ADIC_ARCH) 

all: 

adiC -d gradient -i explosion. init 
g++ -g -o explosion. ad $(AD_INC) $(AD_LIB) \ 
explosion. ad. c explosion. driver. c \ 
-lADIntrinsics-C -laif_grad -lm 

clean: 

rm explosion. ad* 
rm ad_deriv.h 
rm *\~ 

The environment variable ADIC must be set to the directory in which ADIC has been installed. 



C ADIC 1.1 



The following files can be downloaded from 



http : //www-unix .mcs . anl .gov/~naumann/ ad_tools .html 



16 



C.2 explosion. init 



[SOURCE_FILES] 
explosion. c 

[gradient] 
GRAD_MAX=7 

The script file specifies the input hies (only one in this case) and a variety of other parameters 
that can be looked up in The definition of GRAD_MAX, the length of the derivative components, 
is important. Its default value is 5, which would cause trouble in our case where dim=7 > 5. Refer 
to for other ways to set the value of GRAD_MAX. 

C.3 explosion. driver. c 

#include "ad_deriv.h" 
#include <stdio.h> 
#include <iostream> 
#include <stdlib.h> 

extern void ad_explosion(int , DERIV_TYPE* , DERIV_TYPE* , DERIV_TYPE*) ; 

void mainO 
{ 

int i , j ; 

int dim=7; 

int parmax=2 ; 

// independent variables 

DERIV_TYPE *x = new DERIV_TYPE [dim] ; 

DERIV_TYPE *prm = new DERIV_TYPE [parmax] ; 

// independent variables 

DERIV_TYPE *F = new DERIV_TYPE [dim] ; 

// independent variables 

InactiveDouble *jac = new InactiveDouble [dim] ; 

ad_AD_Init (dim) ; 

ad_AD_SetIndepArray (x,dim) ; 
ad_AD_SetIndepDone() ; 

DERIV_val(x[0]) = 1.72; 

DERIV_val(x[l] ) = 3.45; 

DERIV_val(x[2] ) = 4.16; 

DERIV_val(x[3] ) = 4.87; 

DERIV_val(x[4] ) = 4.16; 

DERIV_val(x[5] ) = 3.45; 

DERIV_val(x[6] ) = 1.72; 

DERIV_val (prm [0] ) = 1.3; 

DERIV_val(prm[l]) = 0.245828; 

ad_explosion(dim,x,prm,F) ; 
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for (i=0;i<dim;i++) { 

ad_AD_ExtractGrad(jac,F[i]) ; 
for (j=0; j<dim; j++) 

cout « "f'[" « i+1 « "," « j+1 « "]=" « jac[j] « endl; 

} 

ad_AD_Final() ; 



D ADOL-C 

D.l Modifications of Original Source Code 

#include <adolc.h> 
#include <SPARSE/sparse .h> 
#include <stdio.h> 
#include <iostream> 
#include <stdlib.h> 

extern void explosion_ad(int , adouble*, adouble*, adouble*); 

void mainO 
{ 

int i , j ; 
int dim=7; 
int parmax=2 ; 
int tag = 1; 

// independent variables (passiv) 
double *v = new double [dim+parmax] ; 
// dependent variables (passiv) 
double *Fp = new double [dim] ; 

// independent variables (active) 
adouble *x = new adouble [dim] ; 
adouble *prm = new adouble [parmax] ; 
// dependent variables (active) 
adouble *F = new adouble [dim] ; 



v[0] 




1 


72 




v[l] 




3 


45 




v[2] 




4 


16 




v[3] 




4 


87 




v[4] 




4 


16 




v[5] 




3 


45 




v[6] 




1 


72 




v[7] 




1 


3; 


v[8] 







245828; 
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trace_on(tag) ; 

for(j=0; j<dim;j++) 

x[j] «= v[j] ; 
f or ( j =0 ; j <parmax ; j ++) 

prm[j] «= v[j] ; 

explosion_ad(dim,x,prm,F) ; 

for(j=0; j<dim; j++) 
F[j] »= Fp[j] ; 
trace_of f () ; 

> 

D.2 Computation of Sparsity Pattern 

// Sparsity pattern declarations 
int option [3] ; 

unsigned int** Jsp = new unsigned int* [dim]; 

for(j=0; j<dim; j++) 

Jsp[j] = new unsigned int [dim+parmax] ; 
option [0] =0; // automatic detection for AD mode 
option [1] = 0; // save propagation of bit-pattern 
option [2] =0; //no output 



// Sparsety pattern computation after trace_off 
jac_pat (tag, dim, dim+parmax, v, NULL, NULL, Jsp, option) ; 

D.3 Jacobian Calculation Using Low-Level Routines 

// Calculate Jacobian using forward mode driver 
double** X = myalloc (dim+parmax, dim+parmax) ; 
// Calculate Jacobian using reverse mode driver 
double** U = myalloc (dim, dim) ; 



// Use low level routines to compute Jacobian 
// forward: 

for (j=0; j <dim+parmax ; j++) 
{ 

for (i=0; i<dim+parmax; i++) 

X[j] [i] = 0.0; 
X[j] [j] = 1.0; 
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> 



// first order vector mode forward 

// * ~ = fov_forward 

f ov_f orward (tag , dim , dim+parmax , dim+parmax , v , X , Fp , J) ; 

// reverse 

for (j=0; j<dim; j++) 

{ 

for (i=0;i<dim;i++) 

U[j] [i] = 0.0; 
U[j] [j] = 1.0; 

} 

// prepare reverse sweep with appropriate forward sweep 
zos_f orward (tag, dim, dim+parmax, 1 , v,Fp) ; 

// first order vector mode reverse 

// ~ ~ = fov_reverse 

f ov_rever se (tag , dim , dim+parmax , dim , U , J) ; 



E TAPENADE 

The following files can be downloaded from 

: //www-unix.mcs . anl .gov/~naumann/ ad_tools .html 

explosion. f 
explosion . driver . f 
makefile 

E.l makefile 

ADSTACK = $ (HOME) /ADTOOLS/TAPENADE/tapenadel . 0/stack/ adStack. o 
all: 

tapenade -head expl -vars "x prm" -cl explosion. f 
g77 -o explosion. ad explcl.f explosion . driver . f \ 
$ (ADSTACK) 

clean: 

rm explosion. ad 
rm explcl.f 
rm -fr diffgen 
rm *\~ 

E.2 explosion. driver. f 

program main 
implicit none 
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Example: Explosion Equation 

Driver for computing Jacobian using TAPENADE's 
Reverse Mode 



integer dim, parmax, n 

parameter (dim=7, parmax=2, n=9) 

integer i,j 

independent variables 

double precision x(dim) , prm(parmax) 

derivative components of independent variables 

double precision g_x(dim) 

double precision g_prm (parmax) 

dependent variables 

double precision f (dim) 

derivative components of dependent variables 
double precision g_f (dim) 
the whole Jacobian matrix 
double precision jac(dim,n) 

Initialization of input variables 

x(l) = 1.72 

x(2) = 3.45 

x(3) = 4.16 

x(4) = 4.87 

x(5) = 4.16 

x(6) = 3.45 

x(7) = 1.72 

prm(l) = 1.3 
prm(2) = 0.245828 

Compute Jacobian as sequence of Jacobian transposed 
times vector products 
do 40 i=l,dim 

do 10 j=l,dim 
g_f (j)=0.d0 
g_x(j)=0.d0 

continue 

g_f (i)=l.d0 

g_prm(l)=0.d0 

g_prm(2)=0.d0 

call explcl (dim, parmax, x,g_x,prm,g_prm,f ,g_f ) ; 
do 20 j=l,dim 

jac(i,j)=g_x(j) 
continue 

do 30 j=l, parmax 

j ac ( i , j +dim) =g_prm ( j ) 
continue 
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40 continue 

C print Jacobian 
do 60 i=l,dim 
do 50 j=l,n 

print*, ":£'(", i, j, ") = ", jac(i.j) 

50 continue 
60 continue 

end 
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